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At present there are still several open questions about the origin of the ultra high energy cosmic rays. However, great progress 
in this area has been made in recent years due to the data collected by the present generation of ground based detectors like the 
Pierre Auger Observatory and Telescope Array. In particular, it is believed that the study of the composition of the cosmic rays as 
a function of energy can play a fundamental role for the understanding of the origin of the cosmic rays. 

The observatories belonging to this generation are composed of arrays of surface detectors and fluorescence telescopes. The 
duty cycle of the fluorescence telescopes is ~ 10 % in contrast with the ~ 100 % of the surface detectors. Therefore, the energy 
calibration of the events observed by the surface detectors is performed by using a calibration curve obtained from a set of high 
quality events observed in coincidence by both types of detectors. The advantage of this method is that the reconstructed energy of 
the events observed by the surface detectors becomes almost independent of simulations of the showers because just a small part of 
the reconstructed energy (the missing energy), obtained from the fluorescence telescopes, comes from simulations. However, the 
calibration curve obtained in this way depends on the composition of the cosmic rays, which can introduce biases in composition 
analyses when parameters with a strong dependence on primary energy are considered. In this work we develop an analytical 
method to study these effects. We consider AMIGA (Auger Muons and Infill for the Ground Array), the low energy extension of 
the Pierre Auger Observatory corresponding to the surface detectors, to illustrate the use of the method. In particular, we study 
the biases introduced by an energy calibration dependent on composition on the determination of the mean value of the number of 
muons, at a given distance to the showers axis, which is one of the parameters most sensitive to primary mass and has an almost 
linear dependence with primary energy. 
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1. Introduction 

The cosmic ray energy spectrum extends over more than 
eleven orders of magnitude in energy (from below ~ 10® eV 
to above ~ 10^° eV). It can be approximated by a broken power 
law with four spectral features: the knee at a few 10'^ eV IjJ-lst], 
the ankle at ~ 4 x 10*^ eV iSfil , the cutoff or suppression at 
-3x10'® eVlU^il, and a second knee at - 10'® eV, recently 
reported by the KASCADE-Grande Collaboration ifisll . 

Several experimental techniques are used for the observa¬ 
tion of the cosmic rays, depending on the energy range under 
consideration. In particular, the direct observation of the pri¬ 
mary particles is possible up to energies of the order of - 10'^ 
eV. For larger energies the study of cosmic rays is done by 
observing the atmospheric air showers that they generate as a 
consequence of their interactions with air molecules in the at¬ 
mosphere. There are two classes of ground-based detectors, 
surface detectors and fluorescence telescopes. The surface de¬ 
tectors observe the lateral distribution of the showers by sam¬ 
pling the secondary particles that reach the Earth’s surface, 
whereas the fluorescence telescopes observe the fluorescence 
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and Cherenkov photons generated, during the longitudinal de¬ 
velopment of the showers, as a result of the interaction of the 
secondary charged particles with the air molecules 

Despite great experimental effort done in the last years the 
origin of the cosmic rays is still unknown. The observations 
used to study their origin comprise: the energy spectrum, the 
distribution of the arrival directions, and the composition profile 

aa. 

Certainly, the detailed study of the composition as a func¬ 
tion of energy is of great importance to unveil the origin of 
the cosmic rays at all energies (see Ref. ll^ for a review on 
composition). In particular, it is believed that the composition 
information is crucial to And the transition between the galac¬ 
tic and extragalactic components of the cosmic rays (see for 
instance Ref. and to elucidate the origin of the suppres¬ 
sion at the highest energies 11221] . This feature of the spectrum 
could originate as a result of the propagation of the cosmic rays 
in the intergalactic medium, or by the end of the efficiency of 
the extragalactic sources to accelerate particles at the highest 
energies, or by a combination of both effects. 

At the highest energies (E > 10'^ eV), the composition of 
the cosmic rays is studied by using different observable pa¬ 
rameters obtained from shower measurements which are very 
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sensitive to the primary mass. The parameters most sensi¬ 
tive to primary mass are the atmospheric depth at which the 
maximum development of the showers is reached, which can 
be reconstructed from the fluorescence telescope data, and the 
muon content of the showers or a parameter closely related to 
it, which can be obtained from dedicated muon detectors (see 
for instance Ref. 

The fluorescence light emitted during the shower develop¬ 
ment is proportional to the deposited energy. A fraction of these 
photons is detected by the fluorescence telescopes, making pos¬ 
sible the reconstruction of the longitudinal profiles which yields 
an estimator of the energy of the primary particle. The en¬ 
ergy reconstructed in this way is largely independent of sim¬ 
ulations, just a small correction is done, by using simulations 
of the showers, which corresponds to the so-called invisible en¬ 
ergy 1124 l25n . In contrast, the energy calibration of surface de¬ 
tectors alone has to be done by using detailed simulations of 
the showers and the detectors. The use of simulated showers 
introduces large systematic uncertainties because the hadronic 
interactions at the energies of the cosmic rays are unknown. 
Then, the models used for shower simulation extrapolate low 
energy data taken from accelerator experiments by several or¬ 
ders of magnitude in order to reach the energy of the cosmic 
rays. Note that, at present, the hadronic interaction models are 
being updated in order to reproduce the Large Hadron Collider 
data, which reaches up to cosmic ray energies of the order of 
~ 2 X 10'® eV {Ecn, = 7 TeV) 

The duty cycle of the fluorescence telescopes is ~ 10 % and 
that of the surface detectors is ~ 100 %. Therefore, in order 
to have a large duty cycle and an energy calibration almost 
independent of simulations, the present generation of cosmic 
ray observatories, the Pierre Auger Observatory in the southern 
hemisphere and Telescope Array in the northern hemisphere, 
combine both techniques. The energy scale of the surface de¬ 
tectors is obtained by using a parameter which is in general the 
interpolated signal at a given distance to the shower axis. The 
calibration curve that relates this surface parameter with the pri¬ 
mary energy, reconstructed from the fluorescence telescope, is 
obtained experimentally from a subset of high qu ality events 
observed in coincidence by both types of detectors 111 llll2|] . 

Therefore, if the surface parameter used as an energy estima¬ 
tor depends on primary type, the energy scale depends on the 
composition of the cosmic rays. The use of this energy scale 
in composition analyses introduces biases that can be impor¬ 
tant when the parameters used to infer the primary mass haye a 
strong dependence on primary energy (see Ref. 1341 for a first 
study). 

In this work we study the effects of using an energy scale de¬ 
pendent on composition in mass composition analyses. For that 
purpose we deyelop a dedicated analytical method. We consider 
the AMIGA project ll35[l to illustrate in a simplified but realistic 
way the use of the method. The parameter sensitiye to the pri¬ 
mary mass is the number of muons at 600 m from the shower 
axis 1^ 371, which depends almost linearly with primary en¬ 
ergy. Therefore, composition analyses based on this parameter 
are supposed to be quite affected by the use of a composition 
dependent calibration curye. 


It is worth mentioning that many surface parameters like S * 
S, the risetime of the signals in surface detectors ifssl 1^ . 


the slope of the lateral distribution function Il38l 14011 . the cur- 
yature ratio of the shower front |39, etc. were proposed 
and sometimes used in composition analyses. The composi¬ 
tion analyses that make use of these parameters together with a 
calibration curye dependent on composition are affected by the 
effect studied in this work for the case of the number of muons 
in the context of AMIGA. Each particular case inyolying differ¬ 
ent mass sensitiye parameters and energy calibration methods 
has to be analysed in detail in order to estimate the importance 
of the biases introduced by this practice. 

2. Numerical approach 

2.1. Simulations 

The simulations used in this work are the ones generated 
for the studies done in Ref. S, they correspond to AMIGA. 
AMIGA will consist in a triangular grid of 750 m spacing com¬ 
posed by pairs of detectors, a water-Cherenkoy tank and a 30 
m^ muon counter buried underground. The energy region un¬ 
der consideration goes from lO'^ ® eV up to lO'^ ® eV. 

The atmospheric air showers used in Ref. to pro¬ 
duce the simulated data were generated by using AIRES 1^ | 
with QGSJET-II-03 |^| as the high energy hadronic inter¬ 
action model. The showers were simulated with fixed ener¬ 
gies from log(£’/eV) = 17.6 to log(£’/eV) = 18.5 in steps of 
Alog(£’/eV) = 0.1. In this work just proton and iron primaries 
at 30° zenith angle are considered. Eor each primary energy and 
primary type a set 100 showers were generated. 

The muon counters are segmented in 192 scintillation strips. 
Part of the light generated by a charged particle that goes 
through a giyen strip is collected by a wayelength shifter fi¬ 
bre optic and transported to a multi-anode photomultiplier. The 
electronics has an independent channel for each strip which 
giyes a digital one as output for each muon pulse. Note that 
more than one muon arriying in a time interyal corresponding 
to the width of a muon pulse is counted as one. That is called 
pile-up effect. 

As described in Ref. 1^ . a simplified simulation of the 
muon counters is performed. Muon counters of 100 % of ef¬ 
ficiency buried underground at 2.5 m depth, which corresponds 
to a muon energy threshold of ~ 0.82 GeV, are considered. The 
pile-up effect is also included in the simulations. 

Erom the data taken by the Cherenkoy detectors the lateral 
distribution function of the signal, the impact point, and the 
arriyal direction of the showers are reconstructed. While the 
muon lateral distribution function is reconstructed from the data 
taken by the muon detectors, following the method deyeloped 
in Ref. ll36ll . Note that the reconstruction method includes a 
correction for the pile-up effect. 

The parameter used as energy estimator, considered in this 
work, corresponds to the interpolated signal collected by the 
water-Cherenkoy detectors at 600 m from the shower axis, S , 
and the parameter sensitiye to primary mass corresponds to the 
interpolated number of muons at 600 m from the shower axis. 
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N, a. Note that the signal S is in units of VEM (vertical 
equivalent muon), which corresponds to the average signal de¬ 
posited by a vertical muon that crosses the tank by its center 
14511 . Also note that is the total number of muons in an area 
of 30 m^ for showers at 0 = 30°, therefore, it has no units. 

The AMIGA muon detectors were designed for zenith an¬ 
gles from 0° to 45°. Therefore, we consider here showers of 
30° because it is the median of the zenith angle distribution. 
The parameter for showers of different zenith angles can be 
transformed to 30° by using the corresponding muon attenua¬ 
tion curve H. 

Figure[T]shows versus S for proton and iron primaries for 
three different values of primary energy. 



S [VEM] 


Figure 1: versus S for E = 10'^ ®, lO'*, and eV con'esponding to 

proton and iron primaries of 30° of zenith angle. The high energy hadronic 
interaction model used is QGSJET-II-03. 


The particle density at a given distance to the shower axis for 
the different types of particles of the showers presents shower to 
shower fluctuations. The corresponding distribution functions 
present asymmetric tails. However, the fluctuations introduced 
by the detectors and the effects of the reconstruction methods 
make a Gaussian function a good approximation of these distri¬ 
bution functions. Therefore, the combined distribution function 
of and S , given the primary energy and the primary type, 
can be approximated by a two-dimensional Gaussian function 
which is written as. 


P{N.,S\E,A) = 


1 


In (r[NfA o-[5] Vl - 
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where ((A^), cr[A^]) and ((S ), cr[5]) are the mean value and the 
standard deviation of and S , respectively. The correlation p 


is given by, 

_ cov(A^,5) 

^ cr[A^]cr[5]’ 

where cov(A^, 5) is the covariance between A^ and S. Note that 
(Nfi), {S), cr[A^], cr[5'], and p are functions of primary energy 
(E) and primary type (A). 

The parameters ((A^), cr[A^]) and {{S}, cr[5]) for each pri¬ 
mary energy and primary type are obtained by fitting the corre¬ 
sponding one dimensional distributions with one dimensional 
Gaussian functions. The correlation p is obtained from the 
sample covariance and sample variances corresponding to each 
parameter (see Eq. Q). Figure shows the one dimensional 
Gaussian fits to the proton and iron distributions of Nf^ (top 
panel) and S (bottom panel) for two values of primary energy: 
E - 10*^ ® eV and 10*® ^ eV. It can be seen from the figure 
that the Gaussian fits are a good description of the distribution 
functions. 




Figure 2: One dimensional Gaussian fits to the proton and iron distributions 
of (top panel) and S (bottom panel). The primaiy energies considered are: 
E = 10*^ * eV and 10'*'5 eV. 

In order to obtain an analytical representation of 
P(A^,5|£’,A) the logarithm of the mean value as a func¬ 
tion of the logarithm of primary energy and the standard 
deviation as a function of the logarithm of the primary energy, 
for each primary type and for each parameter (A^ and S ), are 
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fitted. A linear function is used in all cases except the one 
corresponding to cr[A^], for which a better fit is obtained with 
a quadratic function. Also, the logarithm of the correlation as 
a function of the logarithm of primary energy is fitted with a 
linear function for both types of primaries considered. 

The top panel of figure |3] shows the mean value of 5 as a 
function of the logarithm of primary energy for proton and iron 
primaries. Also shown are the fitted points (the error bars are 
included but they are smaller than the markers). The difference 
between these two curves is nearly constant, it increases from 
~9%to~ 11% in the energy range considered. The bottom 
panel of figure [3 shows the relative error of S which is given 
by, e[5] = cr[5]/(5). As expected, it decreases with primary 
energy for both primaries and it is smaller for iron nuclei in the 
whole energy range. 


muons is larger than the one predicted by QGSJET-II-03, then 
S should be more sensitive to the primary mass. Therefore, in 
order to consider more general cases, the difference between the 
mean value of S for iron and proton primaries can be increased 
artificially. For that purpose let us introduce the parameter 6 
(S > 0) such that, 

{S)(E,pr,S) = {l-6)x{S)(E,pr), 

{S){E,fe,S) ^ {l+6)x{S)(E,fe). (3) 

The standard deviation of S for both proton and iron nuclei is 
not modified. Note that 2x6 corresponds to the fraction of 
the average ({S }(E, fe) + {S }(E, pr))l2 added to the difference 
between the mean value of S for iron and proton primaries, as 
can be seen from the following equation. 




Figure 3: Top panel: Mean value of 5 as a function of the logarithm of primary 
energy for proton and iron primaries. Bottom panel: Relative error on the de¬ 
termination of 5 as a function of the logarithm of primary energy for proton 
and iron primaries. 


There are experimental evidences about a deficit on the num¬ 
ber of muons in simulated showers 11471] . As mentioned before, 
the hadronic interactions at the highest energies are unknown. 
As a consequence, models that extrapolate low energy data, ob¬ 
tained in accelerator experiments, to the energy of the cosmic 
rays are used. The signal in a water Cherenkov detector is the 
sum of the electromagnetic (due mainly to electrons, positrons, 
and gammas) and the muonic components. If the number of 


{S )(£, fe, 6) - {S ){E, pr, 6) = A_<5 >(£) + b >(£), (4) 


where A^{S){E) = {S}(E,fe) ± {S}(E,pr). 

The discrimination power of a given mass sensitive param¬ 
eter, q, can be assessed by the commonly used merit factor, 
which is defined as. 


MF(^) = 


kf)fe {^}pr 


(5) 


where Wdx[q]A is the variance of parameter q for primary type 
A. Figure |4] shows the merit factor of A as a function of the 
logarithm of primary energy for three different values of 6: 0, 
0.05, and 0.1. As expected the merit factor increases with pri¬ 
mary energy (see figure|3ll. Also, from the figure it can be seen 
that for 6-0 the merit factor is smaller than one in the whole 
energy range under consideration, indicating that the discrimi¬ 
nation power of S is quite poor. For increasing values of 6 the 
MF increases such that for b = 0.1 it is larger than 1.5 in the 
whole energy range, reaching values close to 3 for energies near 
E - 10**'^ eV. The figure also shows the merit factor of for 
comparison, which is of the order of the one corresponding to 
S forb^O.l. 



Figure 4: Merit factor of 5 as a function of the logarithm of primary energy for 
three different values of the parameter d: 0, 0.05, and 0.1. The merit factor of 
Njj is also shown. 
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2.2. Analysis of the bias 


The mean value of the number of muons at a given distance 
to the shower axis is a parameter commonly used to infer the 
primary mass of the cosmic rays (see for instance 1 48 ' S), due 
to its large sensitivity to the nature of the primary. In this sec¬ 
tion the effects on the determination of the mean value of N^, 
introduced by the use of an energy scale dependent on composi¬ 
tion are studied, considering a realistic physical situation. Also, 


a simplified case is analyzed in Appendix A for which explicit 
analytical expressions of the most relevant quantities discussed 
here can be obtained. 

The calibration curve relates a parameter used as energy es¬ 
timator (the signal S in this case) with the reconstructed en¬ 
ergy. Let us denote the calibration curve as S^^fErec) where C 
indicates that this function depends on the composition of the 
cosmic rays. 

The flux of a given primary can be written as Ja{E) - 
ca{E) J{E), where J{E) is the total flux and ca{E) is the abun¬ 
dance corresponding to a primary of type A. Therefore, 


P{A\E) = ca{E), 

ca{E) J(E) 
P{E\A) = \ j 


P(E) = 


^dEcA{E) J{E) 
JjE) 

£°dEJ(E)’ 


( 6 ) 

(7) 

( 8 ) 


where P{A\E) is the probability to And a nucleus of type A given 
the true energy, P(E\A) is the energy distribution given the pri¬ 
mary type, and P(E) is the energy distribution of all species. 
Note that Yja ca(E) - 1. 

The combined distribution function of the number of muons 
and the reconstructed energy is given by. 


P(N^,Erec\E,A) = P(N^,S^JErec)\E,A) j^(Erec). (9) 

Ot^rec 

Therefore, the distribution function of given the primary 
type and the reconstructed energy is obtained from Eqs. (|7]i and 


P{NAE,ec,A) = 


M(£, 


1 _ r 

rec-) A) Jo 


dE ca(E) J{E) 


xP(N,,S‘rJE 

rec )\E,A), 


( 10 ) 


where, 

M{Erec,A) = 


dN, 

Jo Jo 


dE ca(E) J(E) 


xPiN^,Si,(E,,,)\E,A), 


( 11 ) 


is the normalization of the distribution function. 

A power law energy spectrum is considered in all subsequent 
calculations. Therefore, the total flux can be written as. 


J(E) = C E-\ 


( 12 ) 


2.2.1. Constant composition 

Let us first consider the simplified case in which there are 
just two nuclear species, proton and iron, and that the proton 
abundance Cp is independent of primary energy. Assuming that 
the calibration curve is given by the mean value of the signal 
the following expression is obtained, 

= C/, (SKErecpr) + (1 - Cp) {S){ErecJe). (13) 


Note that the dependence on composition of the calibration 
curve is given explicitly. 

The signal S corresponding to iron nuclei is larger than the 
one for protons (see figure |3]l. Therefore, from Eq. (fOT l it can 
be seen that the reconstructed energy for iron nuclei is larger 
than the true one and for proton primaries is smaller. This bias 
in energy is translated into a bias in composition analyses when 
a parameter with a strong dependence on primary energy, like 
the number of muons at ground, is considered. 

Eigure|5]shows the distribution functions of Np for proton and 
iron primaries for E - Erec - 10*^ eV and for b = 0 (top panel) 
and b = 0.1 (bottom panel). The proton abundance considered 
is Cp - 0.5 and the spectral index is y = 3.27, which corre¬ 
sponds to the experimental value obtained by The Pierre Auger 
Observatory in the energy range under consideration 05111 . It 
can be seen that the distribution functions of proton and iron 
primaries get closer. Note that the distribution function of iron 
nuclei is more affected than the one corresponding to proton 
primaries. This is because the calibration curve tends to move 
the distribution function of protons to the right and the one cor¬ 
responding to iron nuclei to the left, however, the energy spec¬ 
trum tends to move both distributions to the left. Also, from 
the figure it can be seen that the modification of the distribu¬ 
tion functions is more important for increasing values of b, as 
expected. 

The mean value of Np for a given primary type A, as a func¬ 
tion of the reconstructed energy and for a proton abundance 
CAiE) can be calculated from Eq. (fTOl i. 


(NpXErecA) = 


1 


M(Er. 


f^OO 

— dE dNpNpCAiE) 

,A) Jo Jo 

X J(E) P(Np,SijE,,,)IE,A), (14) 
which, by using Eq. ([T]i, takes the following form. 


{Np)(Erec,A) 


(5C,(£,,,)-(5>(£,A))2 

2(t^[S](E,A) 


X 


{Np){E,A)+p{E,A) 


(T[Np\{E,A) 

tr[ 5 ](£,A) 


x[S^^JE,,,)-{S){E,A)) 


(15) 


Note that the correlation introduces a term which is directly 
added to the mean value of Np as a function of the true energy. 


where C is a constant and y is the spectral index. 
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Figure 5: Distribution functions of for proton and iron primaries of E = 
Erec = eV. Top panel: d = 0 and bottom panel: 6 = 0.1. Solid lines 
correspond to the true energy and dashed line to the reconstructed energy. 


Figure 6: Mean value of as a function of the logaiithm of the energy for 
protons, iron nuclei, and a mixture of both primaries such that Cp = 0.5. Top 
panel: d = 0 and bottom panel: 6 = 0.1. Solid lines correspond to the true 
energy and dashed lines to the reconstructed energy. 


Therefore, from Eqs. (|7]i, (|9]l, and (fT4l i it can be demonstrated 
that the mean value of as a function of the reconstructed 
energy, corresponding to a mixture of nuclei, is given by, 

(N^KErec) = Yj^N^){Erec,A) CJA{Erec), (16) 

A 

where ojAiErec) - M{Erec,A)l Yja M{Erec,A). Note that for the 
ideal case in which the reconstructed energy is equal to the true 
energy it easy to show that otAiErec) - CA(Erec)- 

Figure|6]shows the mean value of as a function of the log¬ 
arithm of primary energy for protons, iron nuclei, and a mixture 
of both such that Cp = 0.5 in the whole energy range. Also in 
this case, the solid lines correspond to the true energy (E - E) 
and the dashed lines correspond to the reconstructed energy 
{Erec = E). 

From figure |6] it can be seen that the mean value of Np cor¬ 
responding to each primary is affected by the dependence of 
the reconstructed energy on the proton abundance. In particular 
the mean value corresponding to iron nuclei is underestimated 
and the one corresponding to protons is overestimated. This ef¬ 
fect is quite large for the case of b = 0.1. However when the 
mean value of the mixture is considered there is a cancellation 
that makes the difference between the true values and the ones 
corresponding to the reconstructed energy small. 


In order to quantify the difference between the true value of 
{Np) and the one obtained by using the reconstructed energy let 
us introduce the relative bias (see Appendix A for a simplified 
case in which the bias can be calculated explicitly), which is 
defined as. 


Rh{E) - 


{Np){Erec ^ E) 
{Np){E = E) 


(17) 


Note that Rb is positive for the case in which the true value of 
{Np) is smaller than the one obtained by using the reconstructed 
energy. 

Figure Q shows the relative bias as a function of proton abun¬ 
dance for E - 10'* eV and for b = 0 and b = 0.1. Both curves 
present a maximum between Cp - 0.5 and Cp - 0.7. The rel¬ 
ative bias corresponding to b = 0 takes values between -1.5 % 
and -0.5 %, whereas the relative bias for b = 0.1 takes values 
between -1.9 % and 0.2 %. Therefore, for b = 0.1 the bias is 
extended in a wider range than for b = 0. In the extreme cases, 
Cp - Q and Cp - I, the bias comes only from the convolu¬ 
tion between the spectrum and the energy uncertainty because 
the composition is pure in both cases. From the figure it can 
be seen that for b = 0 the absolute value of the relative bias 
is slightly larger for Cp = 0. This is due to the larger muon 
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content of iron showers, which is more important than the re¬ 
duction of the absolute value of the bias for iron nuclei coming 
from the smaller energy uncertainty (see bottom panel of figure 
O. For the <5 = 0.1 case the absolute value of the relative bias 
for Cp = 0 is smaller than the one corresponding to Cp = 1. This 
is because for 5 > 0 the relative error of S increases for protons 
and decreases for iron nuclei (cr[5] is kept constant). 



J _^^^^^_I_I_^_I_^^^^^^^^_I_^_L 

0 0.2 0.4 0.6 0.8 1 


Figure 7: Relative bias as a function of proton abundance for E = 10** eV and 
for (5 = 0 and (5 = 0.1. 

Figure [8] shows the relative bias on the determination of the 
mean value of Nft as a function of the logarithm of primary 
energy for 6-0 (top panel) and d = 0.1 (bottom panel) and 
for Cp - 0.2, 0.6, and 0.9. The value of the spectral index used 
is the same as before, y - 3.27. For both values of 6 and for 
all values of Cp the absolute value of the relative bias is smaller 
than ~ 2.8 % in the energy range under consideration. This 
value corresponds to the maximum of the module of the relative 
bias as a function of Cp for E - 10*^'® eV. 

As mentioned before, the bias on the determination of the 
mean value of Np depends on the spectral shape of both pri¬ 
maries. Figure |9] shows the relative bias as a function of the 
logarithm of primary energy for y - 0 and for 5 = 0. Compar¬ 
ing this figure with the top panel of figure [8] it can be seen that 
the relative bias has a quite different shape for different values 
of the spectral index y. However, its absolute value is still quite 
small, less than ~ 1.4 % in this case. 

For the case of y = 0 the decrease of the relative bias with pri¬ 
mary energy is dominated by the slower increase of the mean 
value of Np as a function of the reconstructed energy, corre¬ 
sponding to iron primaries, compared with the one correspond¬ 
ing to the true energy. This is due to the fact that as the energy 
increases the merit factor of S also increases (see figurelH, then 
the bias coming from the dependence of the energy scale on Cp 
is more important for increasing values of energy. The mean 
value of Np for iron primaries is more affected by this effect 
because the fluctuations of Np are smaller than the ones corre¬ 
sponding to protons. For the case of y = 3.27 the increase of 
the relative bias with primary energy is dominated by the faster 
increase of the mean value of Np as a function of the recon¬ 
structed energy, corresponding to proton primaries, compared 




Figure 8: Relative bias on the determination of {N^) as a function of the loga¬ 
rithm of the energy. Top panel: d = 0 and bottom panel: ^ = 0.1. Three values 
of proton abundance are considered, Cp = 0.2, 0.6, and 0.9. 



Figure 9: Relative bias on the determination of {Np) as a function of the loga¬ 
rithm of the energy for d = 0 and 7 = 0 . Three values of proton abundance are 
considered, Cp = 0.2, 0.6, and 0.9. 

with the one corresponding to the true energy. In this case the 
fast decrease of the energy spectrum with primary energy tends 
to move the mean value of Np towards smaller values than the 
true ones but this effect is gradually smaller as the energy in¬ 
creases because the fluctuations of S decrease with primary en- 
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ergy (see figure |3ll. As a consequence, the mean value of 
as a function of the reconstructed energy corresponding to iron 
primaries is smaller than the one corresponding to y = 0 but the 
difference with the true value as a function of primary energy 
is almost constant. For proton primaries the mean value of 
is also smaller compared with the one corresponding to y = 0 
also due to the fast decrease of the energy spectrum. As for the 
case of iron nuclei this effect is less important for increasing 
energy making the mean value of to increase faster than the 
true one. 


2.2.2. Varying composition 

Let us consider the case in which the composition profile de¬ 
pends on primary energy. For that purpose the following shape 
for the proton abundance is assumed. 


1 +tanh(fl log(£'/£'o)) 

CpiE) = --- 


( 18 ) 


It represents a transition from iron nuclei at low energies to pro¬ 
tons at high energies. The transition is given at an energy Eq and 
the speed at which this transition takes place is controlled by the 
parameter a. The larger the values of a the faster the transition 
from iron nuclei to protons. 

The top panel of figure [10] shows the mean value of the num¬ 
ber of muons as a function of the logarithm of primary energy 
for Eq - 10'* eV and a — 1. The calibration curve assumed for 
the calculation is given by Eq. (lOT l but in this case the proton 
abundance is a function of energy (given by Eq. (ITSl l). When the 
reconstructed energy is considered (dashed and dotted lines) an 
energy dependent bias appears. Eor both values of 6 considered 
(d = 0 and 6 - 0.1) the transition from iron nuclei to protons 
becomes slower than in the real composition profile. The bot¬ 
tom panel of figure [T0| shows the corresponding relative bias as 
a function of energy for (5 = 0 and (5 = 0.1. It can be seen that 

in both cases the relative bias takes values between-3 % and 

~ 3.2 %. Note that for (5 = 0.1 the relative bias expands over a 
slightly larger region than for 5 = 0 but the difference is small. 


3. Conclusions 

In this work we have studied the importance of a composition 
dependent energy scale on composition analyses. The method 
pursued dwells on a combined distribution function of the mass 
and energy estimator parameters which allows to analytically 
perform all further analyses and their physical interpretation. In 
this paper we have shown the strength of this approach which 
might be applied to different experimental scenarios with ap¬ 
propriate distribution function. This approach allows a clear in¬ 
sight in the impact of different mass composition mixing, which 
is of paramount importance to understand the cosmic ray spec¬ 
trum, in particular in composition changing regions. 

We have applied the method developed to AMIGA in order to 
exemplify these effects in a realistic experimental context. We 
have taken the number of muons and the signal in the water- 
Cherenkov detectors, both evaluated at 600 m from the shower 
axis, as the mass sensitive and energy estimator parameters, re¬ 
spectively. We have found that the distribution functions of the 



log(E/eV) 


Figure 10: Top panel: Mean value of as a function of the logarithm of 
primary energy. Bottom panel: Relative bias as a function of the logarithm 
of primary energy. The dotted and dashed lines correspond to d = 0 and 6 = 
0.1, respectively. The parameters corresponding to the composition profile are: 
Eq = 10** eV and a = l (see Eq. <181 1. 


number of muons for proton and iron primaries can be mod¬ 
ified when an energy calibration dependent on composition is 
used to reconstruct the energy of the events. However, the rela¬ 
tive bias on the determination of the mean value of the number 
of muons is quite small, of the order of a few % in the whole 
energy range under consideration. This is true even for energy 
estimators with merit factors of the order of the one correspond¬ 
ing to the number of muons. We have obtained the same upper 
limit on the relative bias for the two physical situations that we 
have considered: a constant proton abundance as a function of 
primary energy and a smooth transition from iron to proton pri¬ 
maries at £ = 10'* eV. 

It is worth mentioning that the impact of the use of an energy 
scale dependent on composition in composition analyses has to 
be analyzed in detail in each particular case. The reason for 
that is that the effects introduced by this practice depend on the 
parameter sensitive to the primary mass under consideration, on 
the type of detectors used to observe the air showers, and on the 
methods used to reconstruct the shower parameters. 
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Appendix A. A simplified case 

In the simplified physical situation treated here it is assumed 
that the fluctuations of the parameter used to reconstruct the 
primary energy are negligible, i.e. (t[ 5](£’,A) = 0 for A = pr 
and fe. Considering that the energy calibration is given by the 
mean value of the signal, it is possible to obtain the true energy 
E of each primary as a function of the reconstructed energy 
from the following expression, 

<5)(£,A) = Cp{S){E,ec,pr) + - Cp){S){E,ecJe). (A.l) 

Here a binary mixture of protons and iron nuclei and a constant 
proton abundance are assumed. 

The mean value of the number of muons is obtained, in this 
case, following a similar procedure to the one described in sec¬ 
tion |2]2] 


{NpKE 

rec ) = (cp {Np)(E{E rec-i pr),pr) J{E{E 

rec-i pr)) 

dE 

^ i.^rec-> P^) + (1 “ Cp) 

^^rec 

X {Np){E{ErecJe)Je) J{E{E,ec, fe)) 

X -^{E rec^ /e)j X Jcp J{EiE 

rec 7 pr)) 

X -^(E rec-i pr) + {l-Cp) J{E{E 

rec-i fe)) 

Oil<rec 

x-—{ErecJe)\ . (A.2) 

oErec / 

The mean values of parameters Np and S have an almost lin¬ 
ear dependence on primary energy. Then, in order to further 
simplify the calculation let us assume an exact linear depen¬ 
dence on primary energy of the mean value of both parameters. 


By using Eqs. (IA.lllA.21A.31A.4l i and (IA.51 l the relative bias 
on the mean value of the number of muons as a function of the 
reconstructed energy can be written as. 


Rh 


{R^-Rs){Rl-^-\)Cp{\-Cp) 

{cp +Rp{\- Cp)J (cp -H (1 - Cp)) 


where. 


N: 


fe 


Ru = 


Rs = 


N‘" ’ 

c fe 


cPr ■ 


(A.6) 


(A.7) 

(A.8) 


Here a power law energy spectrum of the form J{E) - C E"’’' is 
assumed. 

Prom Eq. (Ia:^ it can be seen that the relative bias is in¬ 
dependent on primary energy. Also, when the composition is 
pure, the cosmic rays are only protons or iron nuclei (cp = 1 
or Cp - 0), the relative bias is zero, as expected. Moreover, the 
bias also disappears when Rs = 1, Rp = or y = 2. Pigure 
lA.l II shows the relative bias as a function of the proton abun¬ 
dance for different values of y, starting from y = 0 up to y = 3.3 
in steps of Ay = 0.1. The values Rs = 1.1 and Rp - 1.5 are 
used to make the plot, which correspond to the ratios between 
the mean values of each parameter (Np and S) for proton and 
iron at £ = 10'* eV. 



Figure A. 11: Relative bias as a function of proton abundance for different val¬ 
ues of the spectral index staiting from y = 0 up to y = 3.3 in steps of Ay = 0.1. 


{Np)(E,A) = (A.3) 

{S){E,A) = (A.4) 

where Eq is a reference energy. Therefore, the mean value of 
Np as a function of the true energy is given by, 

{Np)(E) = (cp + (1 - Cp) . (A.5) 


Prom figure lATTTl it can be seen that the relative bias is an in¬ 
creasing function of y, it goes from negative values at y = 0 to 
positive values for y > 2. This behavior can be understood from 
the fact that for a given value of the reconstructed energy pro¬ 
ton events come from larger values of the true energy but iron 
events come from smaller values, therefore, for a power law 
energy spectrum the iron events have a larger weight, which 
increases with y moving the mean value towards the one cor¬ 
responding to iron nuclei. In this way, after increasing y suffi¬ 
ciently the relative bias becomes positive. 
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For every j + 2, the relative bias has an extreme point placed 
at an intermediate value of the proton abundance. This extreme 
value is a maximum if y > 2 and it is a minimum if y < 2. 
The expression for the proton abundance corresponding to the 
extreme point can be obtained from Eq. (IA.6I) . which is given 

by. 


C 


ext 

P 


R,rV - ^|R,RV 

r.rT' -1 


(A.9) 


which is valid for Rs + I, Rfi ^ Rs, and y 2. As can be seen 
from figure lATTTl varies very slowly with y, in fact it goes 
from ~ 0.54 at y = 0 to ~ 0.58 at y = 3.3. 
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